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Abstract 

Classical reconstruction methods for phase-contrast to¬ 
mography consist of two stages: phase retrieval and to¬ 
mographic reconstruction. A novel algebraic method 
combining the two was suggested by Kostenko et al. 
(Opt. Express, 21, 12185, 2013) and preliminary results 
demonstrating improved reconstruction com- pared to 
a two-stage method given. Using simulated free-space 
propagation experiments with a single sample-detector 
distance, we thoroughly compare the novel method 
with the two-stage method to address limitations of 
the preliminary results. We demonstrate that the novel 
method is substantially more robust towards noise; our 
simulations point to a possible reduction in counting 
times by an order of magnitude. 

1 Introduction 

With the upcoming of high-brilliance X-ray sources, 
new phase-contrast tomography (PCX) methods have 
gained widespread use [2]. Among these is the free- 
space propagation method [7, 21], with the advantage 
that no additional optical elements such as analyzer 
crystals or gratings are required. Compared to conven¬ 
tional absorption contrast, phase contrast may provide 
adequate contrast at lower dose rates, thus allowing 
segmentation of objects comprising two or more mate¬ 
rials with nearly the same electron density (e.g., [26]). 
For a comparison of variants of free-space propagation 
methods in general see [20]. 

First experiments were performed using holo- 
tomography, requiring the combination of measure¬ 
ments from several sample-to-detector distances [8], 
but today highly successful reconstructions are often 
possible using a single distance, e.g., in the area of pa¬ 
leontology [28]. This is experimentally convenient, but 
also remarkable, as the information content is obvi¬ 
ously reduced. In fact, reconstruction with the single¬ 
distance set-up is typically based on the work by Pa- 
ganin et al. [23], assuming proportionality between the 
absorption and the phase shift. Research into the lim¬ 
itations of this so-called duality method will benefit 
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experimental planning. 

The standard approach to PCX reconstruction is a 
two-stage procedure. In the first stage the phase and 
absorption fields are determined for each projection us¬ 
ing a phase retrieval algorithm. In the second stage, a 
classical algorithm is used to compute a reconstruction 
based on the projection fields. 

Recently, Kostenko et al. proposed a combined ap¬ 
proach [17, 18] which - in several of their simulated 
experiments with noise-free data - performs better in 
terms of reconstruction error. However, for simulated 
noisy data the combined method is outperformed by 
the two-stage method. Both methods are tested on 
simulated data with artificial material indices, and 
therefore it is unclear if the combined method performs 
better than the two-stage method for realistic samples, 
for different material types, and for varying noise lev¬ 
els. 

In this paper we provide a careful numerical simula¬ 
tion study of the combined duality method, comparing 
it to the two-stage approach. The use of regulariza¬ 
tion is key to obtaining high-quality reconstruction, 
and we focus on the use of total variation (TV) [25] 
as the regularizer. Many samples in materials science 
and geoscience comprise discrete objects (grains, hbres, 
cracks), and edges naturally lead to high phase con¬ 
trast in the free-space propagation method. Using a 
polycrystal with small density variations as a phantom 
and with simulated noise, our simulations aim to care¬ 
fully compare the reconstruction capabilities of the two 
methods. 

After a review of a classical PCX reconstruction 
method and the combined method, we present our nu¬ 
merical implementation. In simulations with realistic 
material parameters we perform a comprehensive com¬ 
parison of the two methods with respect to different 
material parameters of increasing difficulty and to the 
robustness towards noise. 

Our results show that the combined method pro¬ 
duces improved reconstructions across the range of 
low to high-absorption materials with small absorp¬ 
tion contrast. Furthermore, as the simulated noise 
level is increased, reconstructions from the combined 
method show much greater robustness to the noise. 
This could be of critical importance in practical ap¬ 
plications where noise is always a concern. 


1 



2 Definitions and forward model 


In this section we briefly review the underlying defini¬ 
tions and models. To simplify the presentation and the 
numerical experiments we consider only 2D problems. 

Scalar functions are denoted with italic, e.g., B or 
(j). Vectors are denoted with subscript e.g., Xy or By. 
Matrices are denoted with bold upper case, e.g., A or 
F, and operators with upper case calligraphic letters, 
e.g., A and T. Throughout, || • || denotes the vector 
2-norm. 

2.1 Definitions 

The Fourier transform of a one-dimensional signal / 
is denoted with T: 

/ oo 

da:. (1) 

-OO 

The independent frequency variable is cd and the com¬ 
plex unit is denoted as f = \f^A. The Radon transform 
7^ of a two-dimensional signal f{xi,X 2 ) is defined by: 

[nf]{t,o) = 

f /rtcos(6>) — rsin(^),tsin(^) + r cos(6>)^ dr. (2) 

Here 0 is the angular variable and t is the transla¬ 
tional variable, perpendicular to line integration di¬ 
rection, in other words the coordinate variable on the 
one-dimensional detector. 

A discrete linear inverse problem can be formulated 
for absorption-based computed tomography (CT) by 
discretizing a 2D object into V by V square pixels with 
pixel values stacked into a vector Uy . The Radon trans¬ 
form is discretized using the line-intersection method 
and represented using a system matrix A with elements 
amn of tho path length of X-ray m through pixel n. 
Letting by denote the discrete projection data, the dis¬ 
crete linear inverse problem can be written 

Ally - by. (3) 

For more details see [5]. 

2.2 Free-space propagation model 

Different experimental PCT setups exist, see [2] for an 
overview. In the present work we focus on the free- 
space propagation method. Figure 1, which in some 
sense is the simplest since it does not require analyser 
crystals, gratings or likewise. The set-up is similar to 
the standard CT scanning set-up, with the added re¬ 
quirement that the X-ray source is partly coherent. 

In free-space propagation PCT phase shifts of X-ray 
waves are magnified as the object-to-detector distance 
increases. Based on intensity measurements at the de¬ 
tector, the task is to reconstruct both the absorption 
index and the refractive index decrement of the object. 

For free-space propagation PCT the measured inten¬ 
sity data in the Fresnel region is related to the mate¬ 
rial properties through a non-linear propagation model. 
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Figure 1: 2D sketch of free-space propagation set-up. 
The intensity profile is shown as function of increasing 
detector distance R. In the contact plane R = 0 only 
absorption contrast is visible. In the near-held region, 
of interest here, both absorption and phase contrast 
contribute to the intensity. 



For a 2D object with absorption index S{xi,X 2 ) and re¬ 
fractive index decrement /3(xi,X2), the corresponding 
projections B(t^0) and (j){t^0)^ here called absorption 
and phase shift, can be modelled as line integrals along 
the X-ray propagation [20]: 


Onr 

Be{t)= —mxe), 

(4) 

O'JT 

Mt) = 

(5) 

Here A is the X-ray wavelength and IZ is the Radon 
transform. Based on the absorption and phase shift, 
the measured intensity is modelled as the squared 

absolute value of the convolution of the transmittance 
To and the Fresnel propagator P^\ 

TeO) = exp - Bg{t) + 

, (6) 


(7) 


(8) 


Here R is the object-to-detector distance, | • | the ab¬ 
solute value and ^ the convolution operator. 

Discretization of the domain and hence the object of 
interest into N x N pixels, gives us discrete versions of 
the material parameters (3 and 5. A single index nota¬ 
tion is introduced n = i-\-{N — l)j, for i,j = 1, 2,..., N, 
so n = 1 , 2 ,..., which gives two column vectors, jSy 
and Sy, with element index n. Using the discrete Radon 
transform (3) we dehne discretized versions of Bq and 
( 1 ^ 9 ' 

(9) 

cf>^ = ( 10 ) 

This leads to discrete versions of the transmittance and 

the Fresnel propagator, and from those the intensity 
jR. 

±V . 


Ty = e-Xp{-By (H) 

I^=\Ty*P^\\ (13) 


2 
















Here exp and | • p denote element-wise operations and 
-k denotes the discrete convolution. In practice the dis¬ 
crete convolution is done in Fourier space. 

3 Reconstruction methods 

In this section we review the two methods that we 
compare, namely, a classical two-stage method and a 
combined method presented by Kostenko et al. Both 
methods are equipped with TV-regularization, which 
is described last. 

3.1 Two-stage method 

For absorption CT the measured intensity data can be 
directly related to the material properties, i.e., the at¬ 
tenuation coefficient, by Lambert-Beer’s law [5]. For 
PCT the intensity measurements are related to the 
material properties through the non-linear propagation 
model Eqs. (4)-(8). The standard reconstruction pro¬ 
cess for PCT is a two-stage procedure, consisting of a 
phase retrieval stage and a tomographie reeonstruetion 
stage. 

Different phase retrieval methods have been sug¬ 
gested and investigated in the literature. An intro¬ 
duction to, and a comparison of, some of them can be 
found in [4, 9, 20]. Following the work by Kostenko 
et al. [16, 17, 18] we focus on the contrast transfer 
function (CTF) method, which is based on an assump¬ 
tion of low absorption, and which is derived from the 
expression of the intensity given in (8). Based on a 
Taylor expansion of the transmittance function (6), a 
CTF model that is linear in B and (j) can be derived in 
Fourier space [8]. 

For the case of a single detector distance and the 
so-called duality version of the CTF model with pro¬ 
portionality constant a = —6//3^ we have 

^ [2a sm{7TXR\ujf) — 2cos{7rXR\uj\‘^)]B{uj) 

^Dirac(^)- (1"^) 

Here (^Dirac is the Dirac delta function and uj denotes 
the spatial frequency. For p being the physical size of a 
detector pixel, the sampling distance becomes Fg = 1/p 
and hence uj G [—^, ^]- Introducing discrete fre¬ 
quency values ujm the CTF method can be formulated 
as a discrete linear inverse problem: 

/^ = Wtsd^., (15) 

Wtsd = -2C^ + 2aS^, (16) 

where and are diagonal matrices with mth di- 
agonal elements cos'^^n and sin?/^^, respectively, with 

m = l,2,...,M. (17) 

For the two-stage method based on the CTF duality 
method we use the name TSD. 

The phase-retrieval stage of this method consists of 
solving (15) for By, followed by an inverse Fourier 


transform and then solving (9) for fSy. Since jSy is as¬ 
sumed to be proportional to Sy for this method, we can 
easily find Sy afterwards. 

Following the phase-retrieval stage is a stage where 
a CT reconstruction method is used to compute the 
object of interest. Typically CT reconstructions are 
carried out by using classical methods such as filtered 
backprojection (FBP) [5] or the algebraic reconstruc¬ 
tion technique (ART) [14]. 

3.2 Algebraic combined method 

Methods that combine the phase-retrieval stage and 
the reconstruction stage have been proposed with dif¬ 
ferent aims. In [3] a filtered backprojection type al¬ 
gorithm is derived and tested while in [17] an alge- 
braie combined model is presented and tested, showing 
promising preliminary results. 

In [17] Kostenko et al. also suggested that the stan¬ 
dard phase-retrieval techniques could benefit from us¬ 
ing the redundancy within an entire sinogram rather 
than just being based on the individual projections. 

The term ’algebraic combined’ in [17] refers to a com¬ 
bination of the linear operators. A, describing the dis¬ 
crete Radon transform, F, the discrete Fourier trans¬ 
form, and Wtsd of the linear phase retrieval method, 
i.e., 

Wacd = ^WtsdFA. (18) 

The discrete Fourier transform is introduced because it 
is computationally convenient to formulate the phase 
retrieval model as a matrix multiplication in Fourier 
space. The algebraic combined method based on the 
duality version of the CTF method works on the re¬ 
sulting linear system 

Ip = WacdUv ■ ( 19 ) 

This combined reconstruction method is called the al¬ 
gebraic combined duality method, ACD. 

3.3 Total variation regularization 

In absorption CT, TV-regularization has been shown 

to be advantageous for objects with piecewise constant 
material parameters [1, 12]. TV-regularization pre¬ 
serves edges while smoothing away noise inside homo¬ 
geneous regions. For the discrete linear problem (3), 
we formulate TV-regularization with regularization pa¬ 
rameter a G as: 

2 ;" = argmin i ||Au„ - bjf + aW ||D„u„|| > . (20) 
- [ n=l J 

Here n is the pixel index, is the total number of 
pixels, assuming a square domain, and TlnUy is the 
local finite difference gradient at pixel n. 

Many objects from materials science which would be 
desirable to analyse with PCT have the property that 
they have approximately piecewise constant material 


3 


parameters, e.g., in the form of grains. With this mo¬ 
tivation Konstenko et ah proposed to incorporate TV- 
regularization into both the TSD and ACD methods, 

i.e., in (15) and (19). 

In the present work we also consider TV- 
regularization for both methods and in the remainder 
of the article by TSD and ACD we refer to the TV- 
regularized problems. In a direct assessment of the 
effect of combining linear operators in ACD, and not 
the effect of TV-regularization itself, we find it most ap¬ 
propriate to employ TV-regularization also in TSD. We 
note that one motivation for ACD is precisely the use of 
regularization which through the combination of linear 
operators regularizes the entire reconstruction problem 
including the phase-retrieval step. This is in contrast 
to the TV-regularized two-stage method where only 
the latter reconstruction step is regularized leaving the 
sensitive phase-retrieval stage unregularized. 

4 Our contributions 

In [17] Kostenko et al. compared TV-regularized ACD 
and TSD and demonstrated improvements obtained by 
ACD in terms of root-mean-square error in most of 
their simulation experiments. We find that their pio¬ 
neering results indicate a large potential for ACD, how¬ 
ever we also point to several aspects in which the pro¬ 
vided numerical evidence of reconstruction improve¬ 
ments by ACD may be improved: 

1. The positive results for ACD are for test images 
with one specific choice of material parameters 
that appears to not be motivated from physical 
materials. Thus, it remains open whether as clear 
improvements can be seen for test images with 
physical material parameters. 

2. The positive results for ACD are for noise-free 
data. In fact in a simulation study with noisy 
data, Kostenko et al. find the TSD to be supe¬ 
rior. Only one noise level is considered so it re¬ 
mains unclear whether the combination of phase 
retrieval and reconstruction stages leads to a more 
noise-robust method. 

3. In their implementation of the optimization al¬ 
gorithm to solve the TV-regularized problem, 
Kostenko et al. [17] describe they stop the iter¬ 
ative algorithm when the relative change in the 
objective function value from one iteration to the 
next is smaller than 10“^. This is an intuitive 
choice, however well-known in the field of opti¬ 
mization to be heuristic and not guarantee close¬ 
ness to the solution. This is because the iterative 
algorithm may occasionally take short steps while 
still far from the solution. This means that we can¬ 
not be sure that the shown reconstructions are in¬ 
deed accurate TV-solutions but may be arbitrary 
intermediate images produced by the iterative al¬ 
gorithm. In fact, this problem might affect their 
conclusion that TSD is more robust to noise than 
ACD. 


In the present work, we address all of these three 
problems. Regarding the third problem, we implement 
in our optimization algorithm a stopping criterion that 
does ensure convergence to the TV-regularized recon¬ 
struction, thereby removing any doubt whether the nu¬ 
merical solution returned by the algorithm is in fact the 
sought-after TV-regularized solution. 

We address the first and second problems by pro¬ 
viding two sets of carefully designed simulation exper¬ 
iments. The first set compares TSD and ACD on test 
images with a range of physical material parameters, 
while the second compares TSD and ACD with respect 
to increasing amounts of noise. 

To make the most direct and fair comparison be¬ 
tween TSD and ACD, we employ TV-regularization 
for both and apply the same optimization algorithm 
with the same stopping criterion. Before proceeding to 
the results of the comparisons, we describe in the next 
section our implementation details. 

5 Implementation 

The system matrix A in Eqs. (9) and (10) is large and 
sparse. This means that it can often be stored in the 
memory of a standard modern laptop. For the ACD 
method the dense matrix F makes the combined sys¬ 
tem matrix dense, thus making it infeasible to store 
in memory. We circumvent the problem by a matrix- 
free implementation, in which the applications of the 
forward operator and its conjugate transpose are done 
without explicitly forming the matrices. 

When solving the reconstruction problems Eqs. (9), 
(10) and (19) we impose TV-regularization (20). Solv¬ 
ing such large-scale problems requires efficient algo¬ 
rithms. We chose to implement the Chambolle-Pock 
(CP) algorithm [24, 6] since it was shown to converge 
faster than, e.g., the FISTA method [6] when solving 
problems of the form (3). Moreover the CP algorithm 
can be well suited in the context of CT [27]. 

Our implementation is mainly based on algorithm 4 
in [27], modified by the adaptive parameter approach 
presented in algorithm 2 in [10]. This modified ap¬ 
proach introduces a primal residual and a dual 
residual for iteration k. As mentioned in [10] these 
residuals can also be used to define a stopping criterion, 
since for the CP algorithm we have that 

lim ||pW||2 + ||dW||2 =0. (21) 

/c—)-oo 

We implemented a stopping criterion of the form 

||pW||2 + ||dW||2<r(|bWf + ||dW||2) (22) 

for a user-defined tolerance r. In all of our numerical 
experiments r was set to 10“^. 

The reconstruction stage of TSD involves multipli¬ 
cation with A and its transpose in each iteration. The 
ACD method requires, in each iteration, additional 
FFTs and multiplications with the diagonal matrix W, 
but both operations are much less computationally de¬ 
manding than the multiplications with A and hence 
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Table 1: Absorption index [3 and refractive index decre¬ 
ment S for 40 keV X-rays for the simulated materials. 
From [19]. 


Material 

/5 

(5 

Polycarbonate (C 16 H 14 O 3 ) 

8.43 • 10“^^ 

1.64 • 10"^ 

Carbon (diamond) 

1.90 • 10"^^ 

4.55 • 10"'^ 

Magnesium 

1.15 • 10"“ 

2.22 • 10"'^ 

Aluminium 

2.32 • 10"“ 

3.37 • 10"^ 

Silicon 

2.68 • 10"“ 

3.01 • 10"'^ 

Iron 

6.42 • 10"® 

9.54 • 10"^ 

Copper 

9.96 • 10"® 

1.06 • 10"® 


the computational overhead in an ACD iteration, com¬ 
pared to TSD, is small. 

All implementations and simulations are carried 
out in Matlab, and the implemented code is avail¬ 
able for download at [15]. We also use the func¬ 
tion phantomgallery from the AIR TOOLS package, 
Version 1.3 [11], and the function parbeam from the 
Projector-Pack package. Version 0.2 [13]. 

6 Simulation results 

We compare the TSD and ACD methods across dif¬ 
ferent material parameters and increasing noise lev¬ 
els. The comparisons rely on simulated data, carefully 
modelled to resemble data from a real physical set¬ 
up. Reconstructions are assessed in terms of achievable 
quality, compared to the ground truth and between the 
methods. 

6.1 Experimental set-up 

The simulated experiments are inspired by materials 
science where mappings of structures on a micrometer- 
scale are desired. For polycrystalline materials the 
structures are made up of grains^ and to mimic this we 
use a phantom which resembles grain-structure consist¬ 
ing of three different materials. The phantom, which is 
shown in Figure 2, consists of one background material 
and grains of two different materials, all of which are 
described by indices [3 and S. Indices of the used mate¬ 
rials are listed in Table 1. If nothing else is mentioned, 
polycarbonate is used for the background material. 

Our set-up allows us to simulate free-space propa¬ 
gation PCT experiments and compare reconstruction 
methods. In our simulations we choose specific settings 
realistic for real physical experiments on a laboratory 
X-ray CT scanner. The chosen parameters are pre¬ 
sented in Table 2. To make the simulated data more 
realistic, Poisson distributed noise is used to perturb 
the measured intensity. 

In addition to the settings in Table 2 the duality 
method requires a qualified guess on the proportional¬ 
ity constant a between /3 and 6. This has for all simula¬ 
tions been chosen as the exact proportionality for the 
grain material with the smallest yd, e.g., for the first 
row in Figure 3 a = —1.95 • 10^. Experimental testing 


Table 2: Parameters used in the simulations. 


Parameters 

Settings 

Object: 

2D and 200 x 200 pixels, pixel size 1 

jlVd. 

X-ray 

Energy 40keV. 

source: 

wavelength A = 0.31 A. 

Photons: 

Vo = 10^ photons incident on ob¬ 
ject, average, per pixel, per projec¬ 
tion. 

Distance: 

R = 0.5m. 

Detector: 

Pixel size of 1/im, 572 pixels. 

Projections: 

360 angles 6> G [0°, 180° [ . 



Xi 


Figure 2: 2D phantom with a background material, 
black, and two different grain materials, gray and 
white. 

with different choices of a have shown that the impact 
of changing <j, within ±20% from the exact value, was 
negligible. 

The regularization parameter was chosen empirically 
for each of the simulations in order to achieve the ’best’ 
possible reconstruction. The ’best’ reconstruction is in 
this work measured by two different means: A relative 
error measure 

E = \\u — i^*||/||i4*||, = original, (23) 

and a visual comparison where sharp edges are 
favoured. In the figures with the reconstructions we 
list the specific regularization parameter choices. 

The reconstructions are visualized as images using 
grey-scale color-range [0.9 • min(i^*), 1.1 • max(i^*)]; in¬ 
tensity values outside are truncated to this range. 

6.2 Effects of material properties 

The simulated phantom is varied with materials rang¬ 
ing from low-absorbing material to higher absorbing 
material, i.e., from low yd to higher yd. The two-grain 
materials are chosen such that they have indices nu¬ 
merically close to each other since distinction between 
similar materials is the more challenging case in practi¬ 
cal applications. Increasing the absorption will violate 
the low absorption assumption, which is part of the 
CTF model derivation, so higher absorption is also ex¬ 
pected to increase the difficulty of the reconstruction 
problem. The reconstructions from our simulated ex¬ 
periments with materials of increasing absorption in¬ 
dex are presented in Figure 3. 
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TSD,/?, E=0.41 xio-^° ACD,/?, E=0.38 



Figure 3: Simulations with different materials. First 
row: non-absorbing and non-refracting background, 
grains of polycarbonate and carbon in diamond form 
(highest yd) and regularization parameters arsD = 
0.01, aACD = 120. Second row: grains of sili- 
con(highest yd) and magnesium and regularization pa¬ 
rameters aTSD = 0.1, aACD = 110. Third row: grains 
of silicon (highest yd) and aluminium and regulariza¬ 
tion parameters arsD = 0.12, aACD = 112. Fourth 
row: grains of copper (highest yd) and iron and regu¬ 
larization parameters arsD = 6, aACD = 5. The error 
measure E is defined in (23). 


For these four experiments the ACD method is gen¬ 
erally seen to produce as good or better reconstruc¬ 
tions than the TSD method, in terms of the error mea¬ 
sure E. The TSD results are all visually more blurry 
and with less sharp edges compared to the ACD re¬ 
sults, even though both methods utilize the same TV- 
regularization method. In the process of choosing the 
’best’ reguralized reconstruction from a series of re¬ 
constructions (not shown here), it became clear that 
the TSD reconstructions were corrupted by artifacts 
and/or noise to a higher degree than the ACD recon¬ 
structions. On the four TSD reconstructions in Fig¬ 
ure 3 this is what causes the reconstructions to be more 
blurred, since we gave more emphasis on the regular¬ 
ization term in order to compensate for noise and arti¬ 
facts. We believe that the artifacts are due to the errors 
introduced by the linearization in the phase retrieval 
stage. 

For the experiments with low absorption, in the first 
row of Figure 3, distinction between the different ma¬ 
terials is clear for both methods. The ACD reconstruc¬ 
tion has sharper edges and a lower error measure than 
the TSD reconstruction. 

For the silicon-magnesium and the silicon-aluminium 
experiments, in the second and third row, materials 
with similar chemical structures are seen to be harder 
to distinguish, as expected. The ACD method again 
produces reconstructions with sharper edges and a 
lower error-measure. For the silicon-aluminium recon¬ 
structions in row three, distinguishing between silicon 
and aluminium is difficult for both methods, and alter¬ 
native methods using measurements from two more or 
more distances could improve these results. 

In the experiments with high absorbing materials the 
reconstructions are highly affected by artifacts, such 
that distinction between background, grain, and arti¬ 
fact is difficult. In addition, the edges in the recon¬ 
structions are more blurry. Error measures does not 
tell the same story as the visual inspection, since they 
are relatively low compared to the low-absorption ex¬ 
periments; this is due to the background material being 
correctly reconstructed. 

The ACD method is computationally more demand¬ 
ing than the TSD method, because a larger number of 
iterations is needed to achieve the same solution accu¬ 
racy. In the cases studied here, 1.4 - 6.5 times more 
iterations were needed for the ACD method when using 
the stopping criterion in (22). 

6.3 Effect of noise 

The reconstructions from simulated experiments with 
a decreasing number of recorded photons, and hence 
increasing noise levels, are presented in Figure 4. The 
error measure E is plotted for increasing Nq (cf. Ta¬ 
ble 2) in Figure 5. Using Otsu’s simple thresholding 
segmentation method [22] on the reconstructions, the 
segmentation errors 





xio-^° ACD, d, E=0.26 


xio-' ACD, d, E=0.25 


TSD, d, E=0.3 


TSD, d, E=0.25 


Es = T^misclassified pixels/^^pixels (24) 


are calculated and plotted against Nq in Figure 6. 
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The TSD reconstructions are seen to deteriorate as 
Nq decreases (and the relative noise increases), where 
the grains of the lowest absorbing material closest to 
the object center become indistinguishable from the 
background - cf. the bottom row in Figure 4. The 
edges become more blurry and misclassification of the 
grains is likely to occur. The error measure increases 
drastically to a limit where the reconstructions are un¬ 
reliable. 

ACD reconstructions show much greater robustness 
to the noise: edges remain sharp, materials can be dis¬ 
tinguished, and the error varies slowly with A^o- For 
the problems with higher relative noise (smaller Nq), 
the polycarbonate grains can be visually hard to dis¬ 
tinguish from the background in the chosen gray-scale 
- though numerically the difference is still distinct as 
validated by the low segmentation errors in Figure 6. 



TSD, d, E=0.81 xio-^' ACD, d, E=0.26 



Figure 4: Simulations for increasing data noise, i.e., 
decreasing number of photons Nq. Non-absorbing and 
non-refracting background. Grains of polycorbonate 
and carbon in diamond form (highest /3). First row 
Nq = 5 X 10^, second row Nq = 10^, third row Nq = 
5 X 10^ and fourth row Nq = 10^, photons per pixel. 
Regularization parameter a for the TSD method from 
top to botom were [0.03,0.05,0.08,0.14] and for the 
ACD method [100,140,170,300]. 
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Figure 5: Relative reconstruction error E for increas¬ 
ing number of photons Nq, i.e., decreasing noise in the 
data. 



Figure 6: Segmentation error Es for increasing number 
of photons he., decreasing noise in the data. 

7 Conclusion 

The simplicity of the geometry in the free-space propa¬ 
gation, one sample-detector distance approach makes it 
attractive for many experiments, including those where 
speed of data acquisition or dose is a limitation. Like¬ 
wise, experience has shown that total variation (TV) 
regularization works well for absorption or phase con¬ 
trast tomography on a large class of materials compris¬ 
ing disjunct phases, cracks, pores, etc. 

The outcome of the simulations performed here is 
that the number of photons required to compute a re¬ 
construction of a certain quality can be reduced sub¬ 
stantially and the combined reconstruction method is 
therefore of general interest to the X-ray (and neutron) 
imaging community. We emphasize that the suggested 
combined method is more computationally demanding 
than the classical two-stage method, and hence less 
suitable for on-line or real-time processing. 
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